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^ ■ Abstract 

00 \ We present a numerical study of the time evolution of perturbations of ro- 

(<— ^ \ tating black holes. The solutions are obtained by integrating the Teukolsky 

CN ' equation written as a first-order in time, coupled system of equations, in a 



o 



form that explicitly captures its hyperbolic structure. We address the nu- 
Q^ . merical difficulties of solving the equation in its original form. We follow the 

CJ ! propagation of generic initial data through the burst, quasinormal ringing and 

power-law tail phases. In particular, we calculate the effects due to the rota- 
tion of the black hole on the scattering of incident gravitational wave pulses. 
These results may help explain how the angular momentum of the black hole 
aff'ects the gravitational waves that are generated during the final stages of 
black hole coalescence. 
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I. INTRODUCTION 

Numerical relativity, i.e. the numerical construction of solutions to the Einstein equa- 
tions, is a rapidly advancing field. Reliable mult i- dimensional simulations of astrophysically 
relevant scenarios should become possible in the near future. The systematic exploration 
of the non-linear content of Einstein's theory is of major importance for this accelerated 
progress. In the pursuit of numerically solving non-linear systems, intrinsic accuracy tests, 
such as convergence, help establish the reliability of the solution and its proximity to the 
continuum limit. However, recurrent experiences suggest that tests, extrinsic to the dis- 
cretization process, are also valuable tools for assessing the reliability of complicated nu- 
merical algorithms. One versatile tool, providing a number of such extrinsic tests, is the 
investigation of the linear regime of non-linear systems. In the appropriate limit, numeri- 
cal solutions of non-linear systems should agree with the considerably simpler perturbative 
solutions. This testing approach is bound to have fundamental importance to numerical 
relativity as the field progresses. 

Perturbative methods enjoy renewed attention in their own right. The emphasis now is 
on the application of the standard perturbative methodology to particular realistic systems 
that are thought to be of interest to the emerging field of gravitational wave astronomy. 
Perhaps the most important one among those systems is the collision of inspiraling black 
hole binaries. The key premise of the perturbative approach to black hole collisions is that, 
during the last stages of the coalescence (the close limit), the spacetime can be accurately 
approximated as that of a single perturbed black hole. A particular subclass of such problems 
is the study of black hole collisions as perturbations of Schwarzschild spacetimes; these are 
situations (e.g. head-on collisions) for which the end product of the collision is a slowly or 
non-rotating black hole. Mathematically they are described by the gauge invariant Zerilli 
function |I| . The numerical solution of the Zerilli equation does not pose difficulties and is 
routinely employed in estimating aspects of gravitational wave physics in the context of black 
holes 10]. The generalization of this framework to rotating black holes is physically relevant; 
it is likely that, during the inspiral collision of black holes, the system will not be able to 
radiate all of its angular momentum and will leave behind a single, rotating black hole. 
The work in this paper constitutes an important step towards the goal of generalizing the 
close-limit treatment of black hole collisions to that of perturbations about a single, rotating 
black hole. The general framework to achieve this goal requires: (A) the construction of 
appropriate initial data describing two orbiting black holes, (B) the subtraction of the Kerr 
background from the initial data, to read off the initial perturbations, and (C) the evolution 
of these perturbations. In this paper, we concentrate on the last point, namely the evolution 
of gravitational perturbations of the Kerr spacetimes. 

Previously, we investigated the dynamics of scalar fields in the background of rotating 
black holes [0] (from here on referred to as Paper I). We considered both slowly and rapidly 
rotating black holes. For slowly rotating black holes the background geometry can be treated 
as a perturbed Schwarzschild spacetime, with the angular momentum per unit mass playing 
the role of a perturbative parameter. The study was centered on the late-time dynamics of 
the scalar field and showed that, for rotating black holes of arbitrary angular momentum, the 
late-time dynamics is dominated by the lowest allowed multipole (/ = 0). Here we undertake 
the development of a numerical scheme for the study of gravitational perturbations of Kerr 



black holes in the time domain. The main objective of this work is to follow the propagation 
of generic initial data through the burst, quasinormal ringing and power-law tail phases 
of the evolution. In particular, we are interested in investigating the effects of black hole 
rotation on the scattering of incident gravitational wave pulses. A characterization of such 
effects provides an indication of the role that rotation plays on signals produced during the 
final stages of black hole coalescence. 

II. PERTURBING KERR BLACK HOLES 



At first instance, a direct derivation of the equations governing the perturbations of 
Kerr spacetimes is to consider perturbations of the metric. This path, however, leads to 
gauge-dependent formulations. A theoretically attractive alternative is to examine curvature 
perturbations. Using the Newman- Penrose formalism, Teukolsky [|,Q derived a master 
equation governing not only gravitational perturbations (spin weight s = ±2) but scalar, 
two-component neutrino and electromagnetic fields as well. In Boyer-Lindquist coordinates 
and with the use of the Kinnersley null tetrad M, this master evolution equation reads 
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where M is the mass of the black hole, a its angular momentum per unit mass, and A = 
r^ — 2Mr + a?. For gravitational perturbations, the function \l/ is given in terms of the 
Weyl tensor tetrad components ipo and ifj^. That is, ^ = ■i/'o for s = +2 and \l/ = p~^iIja for 
s = —2, with p = — l/(r — mcos6'). 

The Teukolsky equation (|l|) reduces to the Bardeen-Press equation [0 in the limiting 
case of non-rotating black holes. For the case s = 0, it yields the equation for a scalar wave 
propagating in a Kerr background, a system which we studied numerically in Paper I. 

A result of great importance for perturbation theory was the discovery by Teukolsky 
that, when viewed in the frequency domain, Eq. (|lD is separable in the r and 9 coordinates 
(separation of the azimuthal angular dependence is always possible). To our knowledge, 
most of the work on the dynamics of perturbations of Kerr spacetimes has been performed 
in the frequency domain (or under the assumption of a harmonic time dependence). This has 
certainly been the case for studying quasinormal mode frequencies |]8|-[TI|, wave scattering 
|l^ , and the motion of test particles in the Kerr background |[T3HT5| . 

Here we are interested in the time integration of physical initial data, possibly from the 
inspiral collision of binary black holes. The data may be analytic approximate solutions to 
the linearized constraints, or numerical solutions of the non-linear constraints. 

In principle, one can Fourier transform the initial data and perform the evolution of such 
data for each frequency. Afterwards, the data are transformed back to the time domain if 
needed. From the computational point of view, however, this approach is far more expensive 
than it is to solve the problem in the time domain. For several reasons the number of 



frequencies that one has to consider is orders of magnitude higher than the number of angular 
components required to resolve the 6'-direction: While the angular directions are bounded, 
the time direction is only bounded from below by the initial data surface. Focussing on the 
investigation of quasinormal modes (QNM) and power-law tails, one would first of all expect 
that one needs quite a fine resolution near the u = point in order to be able to correctly 
resolve the tails. Similarly, the resolution of the QNM would be very sensitive to the spacing 
in frequencies. With this in mind, we have chosen to solve the Teukolsky equation in t, r, 
and 9 coordinates. 

The resulting 2+1 evolution equation is a hyperbolic, linear equation which is quite 
amenable to numerical treatment, provided suitable coordinates, variables and boundary 
conditions are chosen. The major numerical difficulty in solving the Teukolsky equation in 
the time domain arises from the term linear in s, involving the first-order time derivative. 
Depending on the relative sign between the coefficients of dt"^ and du"^, one may view dt"^ 
either as a damping term (when the signs of both coefficients agree) or an anti-damping 
term (otherwise). Without the factor 2s, the real part of the coefficient of dt"^ reads 

C{r) = -^ '- - r. (2) 

In the physically allowed range [r_|_, oo), where r+ = M + ^ M"^ — o? represents the event 
horizon, the function C is monotonic in r, with limr^r.+ C = oo and limr__>ooC = —oo. 
Therefore, there exists a point r^ such that C{r < Tc) > and C{r > Tc) < 0. On the 
other hand, the coefficient of dtt^ is always negative. Consequently, for s > 0, the term 
dt^ acts as a (anti-) damping term for (r < r^r > r^.. The situation reverses for s < 0. 
Of course, there is no real damping or anti-damping, in terms of energy balance. The 
precise polynomial dependence of the coefficient (s) is affecting, among other things, the 
required asymptotic fall-off behavior of ingoing and outgoing radiation. However, the above 
consideration suggests that, in numerical work, the presence of such derivative terms could 
lead to exponentially growing modes. Those unphysical modes are not part of the family of 
solutions to the Teukolsky equation. Indeed, in our first attempts of solving the Teukolsky 
equation, such growing modes were present in our simulations. These modes were clearly 
of numerical origin. The suppression of these instabilities proved to be an interesting and 
challenging exercise in the construction of numerical algorithms. 

The two key factors in successfully solving the Teukolsky equations were: first, to rewrite 
equation (|^) in a form that explicitly exhibits the radial characteristic directions, and second, 
to carefully select the evolution field and its asymptotic behavior. 

A successful numerical evolution of the Teukolsky equation was achieved by considering 
the s = —2 case. This restriction does not diminish in any way the range of initial data 
that can be evolved. In retrospect, the choice of ip4 as the evolution field is natural since 
the radiation content at infinity is most succinctly described in terms of that field. Yet 
in our investigation this choice was dictated by the need to obtain numerically acceptable 
asymptotic limits near the horizon. In order to reconstruct the perturbation of the metric 
coefficients from the solution of the Teukolsky equation, however, one would have to evolve 
the equation for s = 2 as well. In our approach, i.e. evolving tp^, it is possible to recover the 
metric perturbations only at infinity. In some sense this restricts the information that we get 
from our calculations. Local information about the perturbed spacetime (about geodesies 
etc) is not available. 



For a given spin weight s, outgoing waves correspond to solutions to Eq. (||) with the 
hmiting behaviour |^ 

hm 1^,1 ~ l/r^'+^ (3) 

r*^+oo 

hm l^^l ~ 1 (4) 
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at spatial infinity and the event horizon, respectively. Meanwhile, ingoing waves behave as 

lim I^J ~ 1/r (5) 

lim 1^,1 ~ A-^ (6) 

Here we have used the Kerr tortoise coordinate r*, which is defined by 

dr* = —^dr . (7) 

The convenient properties of the s = —2 choice can be verified by looking at the asymp- 
totic form of propagating waves near the horizon (r* -^ — oo). For s = —2, the solutions 
are bounded for any direction of propagation; in contrast, for s = 2 the ingoing solution 
diverges as the horizon is approached (A -^ 0). The asymptotic behavior for s = —2 at 
r* — > +00 can be subsequently fixed by requiring that outgoing waves have asymptotically 
bound amplitude also in this limit. For the Teukolsky function \&_2, this is achieved through 
a rescaling by an appropriate function of r. A convenient and simple choice is r^, a factor 
that is regular at the horizon. 

Regarding the choice of spatial coordinates, we use the Kerr tortoise coordinate r* defined 
by Eq. (0). As azimuthal coordinate, we use the Kerr coordinate instead of the Boyer- 
Lindquist coordinate (p. The coordinate transformation to (f) is defined by 

d^ = d(t)+ ^dr . (8) 

The azimuthal coordinate was previously used for the scalar field evolution in order to 
ameliorate coordinate pathologies near the horizon. A discussion of those pathologies and 
their precise manifestation in the slow rotation limit is given in Paper I. 

III. SOLVING THE TEUKOLSKY EQUATION 

We can now introduce the following Ansatz for the solution to the Teukolsky equation: 

^{t,r*,e,4))=e'"'^r^^t,r\e) . (9) 

It follows that the Teukolsky equation for $(t, r*,6') has a structure similar to that of Eq. 
(0), and therefore the same analysis of the damping or anti-damping nature of first-order 
time derivative terms applies. After a series of unsuccessful numerical experiments with this 
second-order in time and space equation for $, we found that numerical instabilities were 
suppressed by introducing an auxiliary field 11 that converts the Teukolsky equation to a 
coupled set of first-order equations in space and time. This can be accomplished by defining 
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The resulting form of the Teukolsky equation has the following structure: 

dtu + Mdr*u + Lu + Au = 0, 
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where u = {^r, $/, Hr, Hj} is the solution vector with (1) and (R) labeling the imaginary 
and real parts, respectively. The coefficient matrices M and A are given by: 
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respectively. Finally, X is a matrix operator that contains all the angular derivatives and 
has the following non-vanishing elements: 
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The coefficients in the above matrices are given explicitly in Appendix A. 

The first-order system given by Eq. ( [TT| ) is hyperbolic in the radial direction since the 
eigenvalues of M are real and there is a complete set of linearly independent eigenvec- 
tors. Diagonalization of the matrix M and construction of evolution schemes based on the 
eigenfields generated stable evolutions. However, it turned out that this last step was not 
necessary. Stable evolutions were also achieved using a modified Lax-Wendroff method (see, 
for example, |jl6[) applied to Eq. (|ll]) when this equation was rewritten as: 



dtU + Ddr*u = S 



(15) 



with D = diag(6, b, -b, -b) and S = -{M - D)dr*u - Lu - Au. 

Equation (|T^) is discretized on a uniform two-dimensional polar grid. Typically we used 
a computational domain of size -lOOM < r* < 500M and < 6^^ < vr with < i < 8000 
and < j < 32. The Lax-Wendroff updating of field values it" given on a timestep t„ 
consisted of two steps. In the first step, intermediate field values n" -, /(, are obtained from 
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where we have omitted the angular indices j for clarity. The fields are j-centered in the 
angular direction with angular derivatives approximated by the appropriate second-order 
centered expressions. Radial derivatives in S'^_^_i/2 are approximated by centered difference 
expressions of field values at i and i + 1. Similarly, algebraic terms in -D"+i/2 and -S'^^.^^m 
are obtained by averaging over the values at i and i + 1. The second and final step is a 
Staggered leapfrog step: 



ur+i = u?-6t 
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Here the centered difference expressions and averages for S^ and D" are taken from 
values at u^_^_-^^ L and u^_i/2 ■ Examination of radial and angular characteristic directions 
at arbitrary locations shows that the appropriate CFL condition is 6t < max{ 5r*, 5 56'}, 
where 6t is the evolution timestep. 

Boundary conditions are required at the horizon, at the rotation axis of the black hole 
and at radial infinity. The boundary conditions near the horizon are $ = 11 = 0, as follows 
immediately from the asymptotic behavior of ingoing waves (see Eq.(^). The distinction 
of ingoing and outgoing waves near the Kerr event horizon requires a careful definition 
because of the rotational dragging of inertial frames [Q. However, in practice it clearly 
does not matter if we set both propagation modes to zero. At the axis, only a condition 
on $ is required. Depending on the behavior of the angular eigenfunction belonging to the 
azimuthal number m near the axis, one imposes either $ = or dg^ = |T2[. Finally, 



the boundary conditions at infinity are much harder to impose when using a finite radial 
grid. Approximate schemes allow the transmission of waves across the boundary, but the 
implied truncation of the equations coefficients interferes strongly with physical features 
of the evolution such as quasinormal ringing and tails. A clean solution can be achieved 
via "Cauchy-characteristic matching" (see p!7| , p!8| and references therein). However, this 
approach will not be pursued here. For the problems of immediate interest the computational 
domain can be made sufficiently large that errors generated at the outer boundary will not 
affect the results. Consequently, we set $ and H arbitrarily to zero at the outer boundary, 
and our results are computed only from information inside the maximal Cauchy development 
of the initial data surface. 

The stability of the code was verified with long-time evolutions, of the order of lOOOM. 
Its accuracy was tested using standard convergence tests. The code was found to be second 
order convergent for evolution times t < 50M . The convergence rate degrades as the total 
evolution time is increased, but is consistently above 1.3. A note on the nature of the un- 
stable modes discussed previously is of some relevance here. Linear systems with constant 
coefficients may exhibit numerical instabilities if the eigenvalues of the update matrix have 
modulus larger than unity. This phenomenon is easily analyzed with the so-called Von Neu- 
mann analysis of the difference equation. A similar analysis can be used in equations with 
variable coefficients, via local stability analysis, but the method becomes far less conclusive 
in that case. For the Teukolsky equation, a number of simple, locally stable, discretizations 
turned out to be unstable for late times. The instabilities depended on the numerical dis- 
cretization lengths, with increased resolution generally leading to slower growth rates for the 
unphysical modes. Yet only impractically high resolutions would suppress the instabilities 
for the long evolution times required for the applications discussed here. 



IV. BURSTS, QUASINORMAL MODES AND POWER-LAW TAILS 

The initial data for gravitational perturbations of rotating, asymptotically flat spacetimes 
can be broadly divided into two classes: Waves coming in from infinity and then scattering 
off the black hole background, and waves emanating from the black hole as a result of an 
external excitation of the black hole. The latter case is perhaps the most interesting since 
it includes the description of the final stages of the binary black hole coalescence within the 
close-limit approximation. 

In general, the evolution of both types of initial data consists of three stages, as seen 
by an observer located away from the hole. During the first stage, the observed signal 
depends on the structure of the initial pulse and its reflection by the curvature potential 
(burst phase). This phase is followed by the exponentially decaying quasinormal ringing of 
the black hole (quasinormal phase). In the last stage, the wave slowly dies off as a power-law 
tail (tail phase). The precise manifestation of the last two phases is dictated by the subtle 
interference of the respective amplitudes. The complex frequencies of the QNMs for various 
values of a are well known ISHTTl and we can use them as benchmarks of the present code's 
accuracy. Similarly, the late-time tail phenomenon is mathematically well understood in 



spherically symmetric backgrounds, where tails have been predicted and verified |ll9| , |20 
For rotating black holes, the background is no longer spherically symmetric, but numerical 
results for scalar waves in Kerr backgrounds suggest that power-law tails will also be present 
for spin-2 fields 0]. Furthermore, it has recently been argued that these tails should take 
the same form as for non-rotating black holes pT| . 

We have used the numerical algorithm described in the previous section to obtain the 
time evolution of generic perturbations impinging on the rotating black hole. For all the 
results presented here, we have used ingoing initial data with compact support for $/j and 
have set $/ = 0. All the evolutions clearly showed the exponentially damped oscillations of 
the QNMs and the subsequent late-time power-law tail. We will now discuss each of these 
phases in somewhat more detail. 

A. Power-law tails 

In Paper I we showed that the late time evolution of scalar fields in the background of 
rotating black holes is qualitatively similar to the non-rotating case. Here we extend this 
result to the physically more interesting evolution of a spin-2 field. 

Starting with initial data oc -2^2™ (where s^i^ denotes a spin-weighted spherical har- 
monic), we studied the late time regime for a/M = 0, 0.25, 0.5, 0.75, 0.9 and m = 0, ±1, ±2. 
In all cases, we found a power-law tail behavior |$| oc t^^. For a nonrotating black hole 



it is known that the power-law exponent should be yU = 2/ -|- 3 pO| , |2T| . Indeed, for a = 
we could reproduce the theoretical value /i = 7 with an accuracy of 10%. Moreover, our 
calculations show that the exponents governing the behavior for a 7^ remain similar to 
the Schwarzschild value. These results are summarized in Table |. We conclude that the 
power-law tail behavior is basically determined by the dominant asymptotic form of the 
potential, namely its "Schwarzschild" part. That this would be the case has recently been 
suggested by analytic arguments [^. Analytic results (for non-rotating black holes) can 
also explain the main discrepancy between our numerically determined power-law exponents 

8 



and the theoretical values. The field will be governed by a pure power law t (^'+3) only at 
very late times. In a somewhat earlier regime "higher order" terms will also be significant 



2T[. The presence of such higher order terms tend to increase the value of a numerically 
extracted power-law exponent. 

The results discussed above correspond to initial data that are dominated by the lowest 
allowed multipole / for the particular value of m that is being considered (/ > Iml). When 
this is not the case, we find that different multipoles become mixed in the evolution. To 
investigate this, we started from an initial angular distribution given by -2^3™, and then 
computed the power-law exponent for different values of m. For m = 3 the initial pulse 
used for the two-dimensional evolution corresponds to the lowest possible value of / and 
the late time power-law behavior is consequently dictated by /x = 10.05. The situation 
is different for the case m = 0, where we find that the power- law exponent is given by 
/i = 7.72. That is, the late-time evolution is dominated by the quadrupole. In contrast, the 
corresponding Schwarzschild evolution exhibits no such mixing of multipoles, and evolutions 
for I = 3,m = 0,3 agree well with the theoretical power-law tail exponent /x = 9. A result 
similar to this was discussed in Paper I. 

This result is, however, not too surprising. The mixing of multipoles is due to the rota- 
tional dragging of inertial frames. Basically, there are two reasons why different multipoles 
will be present in the evolution. The first one is "imperfect" initial data: Here we have 
expressed the initial data in terms of the spin-weighted spherical harmonics s^i"^. The 
appropriate angular functions that are used when the r and 6 components of Eq. (|l]) are 
separated are the spin-weighted spheroidal harmonics sSp{9;auj) 0). Hence, because of 
the frequency (time) dependence in the sSp{0] au), it is difficult to generate initial data for 
the Kerr problem that represent a pure multipole / for a specified value of m. Furthermore, 
the rotation of the black hole will lead to an active coupling between different multipoles 
(for a discussion of the analogous situation for rotating stars, see [^). So even if we could 
initiate the evolution with a pure multipole, other multipoles would be generated and may 
play a role at later times. 

The mixing of multipoles can have an interesting result. Suppose that the initial data is 
dominated by a certain multipole / but that a relatively small part of / — 1 is also present. 
Each of these multipoles will give rise to power-law tails, but the amplitude of the latter 
will be much smaller than the first. In this situation the late time field should contain both 
a tail of form ^-(2^+3) g^^^j another term that behaves as t~(^'+^\ Since the first term dies 
off faster than the second it seems unavoidable, no matter how small the amplitude of the 
second tail term is, that the very late time evolution of the field will be dominated by the 
/ — 1 multipole. The corresponding time evolution will simply exhibit a late-time change in 
the angular behaviour. An example of this was discussed in Paper I. 

B. Quasinormal modes 

To further test the performance of the Teukolsky code, we performed a series of simula- 
tions for different values of the angular momentum parameter a and the azimuthal number 
m. Then we focussed on the later part of the ringing sequence, which is dominated by the 
slowest damped mode of oscillation, and read off the oscillation frequency u via a standard 
Fourier transform of the time signal. Table |T| summarizes results obtained in this way. 



A comparison of our values for the real part of each QNM frequency with those given 



by Leaver and Kokkotas |11| for m = 0, a/M = 0, 0.5, 0.9 shows an agreement to better 
than 0.3%. For m = —1, our results agree with the values of Detweiler and Leaver 
(different sign convention) to within 0.5%. Our results for the imaginary parts are accurate 
to 1%. Although higher multipole evolutions are readily obtained, accurate estimates of the 
corresponding QNM frequencies require increasingly higher angular resolution. Hence, we 
have only estimated QNM frequencies for the quadrupole modes. But it should be pointed 
out that, as was the case for the power-law tails, QNMs corresponding to other multipoles 
are present in all signals. 

Interestingly, we find that the numerically extracted QNM frequencies for non-zero m 
do not depend on the sign of m, i.e. we get the same values for the QNM frequencies from 
evolutions for e.g. m = ±1. At first sight this is surprising, since it seems to contradict 
well established results. According to for example Detweiler |P , the imaginary parts of the 
m = +1 and m = —1 modes should become quite different as a increases. 

Fortunately, the answer is simple: The frequencies of both the m and the —m QNMs 
are present in a typical evolution. This also explains the existence of two distinct regimes 
of the quasinormal ringing that can be seen in Fig. Q. In this figure we show the field as 
seen by an observer at r* = 20M, 6 = 7r/4 for a = 0.99M, M = 0.5, and m = 2. In the 
early phase of QNM ringing (120M < t < 160M), the decay of the field is approximately 
given by |<I>| oc e^°'^^^*, whereas the behavior for late times is governed by a slower decaying 
mode, |$| oc e~°'°^°^*. These values are in good agreement with the theoretical values for 
/ = 2 and m = ±2, as can be seen by comparison with Fig. 1(a) in p|. 

One question remains. Why should we expect both these modes to be present in the 
signal? The answer can be found in the symmetries of the problem. In the Schwarzschild case 
we know that there will be modes in the first and second quadrant of the complex ci;-plane 
(assuming that the modes have a time-dependence e*'^*). All these modes will contribute to 
the signal, but since they occur as complex conjugate pairs uji and —ujf (where the asterisk 
represents complex conjugation) the contribution from one set of modes can be expressed in 
terms of the other. Hence, the QNM part of the signal can be evaluated using only modes 
in the first quadrant (for a study of the Schwarzschild problem, see ||21|| ). 

In the case of Kerr black holes the situation is more complicated. The QNMs no longer 
occur as complex conjugate pairs. Let us suppose that we work in the frequency domain, and 
that we have a solution "^imiuj, r, 6)e^"^'^ to the Teukolsky equation for given / and m. Then 
it is easy to show (using the separated equations, see ^) that a solution to the equation for 
—m will follow as 

^i-U^, r, 9) = l^iU-^*, r,n- 9)]* . (18) 

This means that if u; is a QNM corresponding to / and m then the complex conjugate —uj* 
will be a QNM for / and —m. This symmetry is nicely illustrated in Fig. 3 in 0. We can 
use this symmetry to express the contribution from the Kerr QNMs in the second quadrant 
in terms of modes in the first quadrant (we will describe this procedure in detail elsewhere). 
If we represent modes in the first quadrant by uimn, we can schematically write the QNM 
contribution to the signal as 

^T"" = E {Fi^lmn) + G(-<_)} . (19) 

In 
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This explains the presence of both branches of QNMs in our evolutions. 

C. Wave propagation 

Finally, an interesting direct application of our code is the identification of effects of the 
background rotation on the propagation of signals. The idea would be to use the same initial 
data for black holes with different rotation rates and see to what extent the rotation of the 
black hole can be inferred from the emerging signals. A stumbling block for work in this 
direction, though, is the difficulty to prescribe initial data that represent the "same" initial 
physical perturbation for any value of the black hole angular momentum. This prescription 
is possible for data coming in from infinity, since the effects of rotation are negligible for large 
enough radii. However, it is not clear how to set up initial data (e.g. position of the pulse) 
in the vicinity of the black hole in such a way that comparison of evolutions for different 
angular momenta makes sense. For this reason, we addressed only the effect of rotation on 
pulses that were initially far way from the hole. 

In Fig. we show a series of snapshots in r* of the evolution of a pulse for m = 2. At 
t = 0, the initial data consists of a bell-shaped pulse in r* centered at 75M and propagating 
inwards. Then the evolution of $/j is shown at constant 9 = 7i/2 and < t < 200M for 
four different values of the angular momentum of the black hole: a/M = 0, 0.5, 0.9, 0.99. 
One immediately sees that an increase in the angular momentum of the black hole has two 
effects on the scattered pulse. Both the amplitude and wavelength of the emerging signal 
change. By t = 200M the wavelength in the trailing part of the signal for a/M = 0.99 is 
approximately half as long as for the non-rotating case. This agrees well with the anticipated 
increase in the dominating QNM frequency ||^. Furthermore, at t = 200M the amplitudes 
for a/M = and a/M = 0.99 also differ by ~ 25%. The amplitude of the emerging signal 
generally decreases as we increase a. 

The angular behaviour for the scattering scenario also shows interesting features. An 
example can be seen in Fig. ^ where we show the evolution for a/M = 0.99. The initial 
data is given by ^r cx sin^ 9, centered in r* at 75M and propagating inwards, but after some 
time there is a clear transition of the angular distribution to -2^2^ oc cos'^9/2. This is the 
dominant angular behavior during the later evolution. That is, the quadrupole dominates 
the scattered wave at later times. This situation is similar to the one we discussed for 
power-law tails. 

V. CONCLUSIONS 

We have studied the time evolution of perturbations of rotating black holes by numer- 
ically integrating the Teukolsky equation written as a system of equations of first-order in 
time and space. In particular, we investigated the role of the black holes angular momentum 
on the dynamics during the burst, quasinormal ringing and power-law tail phases. We have 
reproduced values for the fundamental quasinormal mode frequencies with an accuracy of 
better than 1%. As for the late-time regime, we found power-law tail behavior analogous 
to the scalar field case that we had studied previously. These results are in good agreement 
with what one would expect and establish the reliability and accuracy of our code. We 



11 



thus have access to a useful numerical laboratory that can be used to further illuminate the 
detailed physics of rotating black holes. 

In the future we plan to use this laboratory to investigate several interesting issues. We 
would, for example, like to establish what effect the so-called super-radiance will have on the 
evolution of initial data. From studies of scattering of monochromatic waves from rotating 



black holes (see, for example |T^) it is know that frequencies such that ujM < am/2r+ 
will correspond to a reflection coefficient whose magnitude is larger than unity. That is, 
these frequencies are "super-radiant" . This is the wave-analogue to the well-known Penrose 
process. In the results presented in this paper the effect of super-radiance was not obvious, 
but it is plausible that it can be unveiled in a more detailed study. For example, it seems 
that super-radiance should lead to different results for co- and counter-rotating waves in an 
evolution. That is, if super-radiance is relevant one should see difference between evolutions 
for |r7i| and — |m|. 

Another interesting issue concerns the excitation of quasinormal modes. It is known that 
the damping of the —\m\ modes will vanish as a increases 0. This would potentially make 
the modes of a rapidly rotating black hole ideal for gravitational wave detection. However, 



in an approximate study Ferrari and Mashhoon [^ have shown that the amplitude of these 
extremely slowly damped modes will be negligible. In the extreme Kerr limit no energy 
is expected to be radiated through these modes. This result has proved to be very hard 
to verify analytically, but it is testable with the present code. We should thus be able 
to establish whether the extremely slowly damped quasinormal modes are of astrophysical 
relevance or not. 

Another outstanding problem in black-hole physics concerns the dynamical stability of 



the Kerr black hole to perturbations. So far, only mode-stability has been established [24 



To prove complete stability one must ensure that all relevant quantities remain point-wise 



bounded during an evolution, cf. [|25[. It seems plausible that our Teukolsky code can prove 
useful also in this context. 

Perhaps the most interesting future application of our numerical Teukolsky code will be 
to evolve initial data qualitative similar to the late stages of a binary black hole coalescence. 
The main problem in approximating black hole collisions with perturbations about Kerr 
spacetimes is the construction of initial data. Work so far on the close limit approxima- 
tion to construct initial data has heavily depended on having a conformally fiat background 
geometry. On the other hand, constant time hypersurfaces in Boyer-Lindquist or Kerr coor- 
dinates are not conformally fiat. Hence disentangling these effects is not as straightforward 
as for nonrotating black holes. Possible ways of circumventing these obstacles are presently 
being investigated [^. 
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APPENDIX A: 



In the following, we give the explicit form of the coefficients in the matrices M, A, and 
L in Eq. (|Tl|). With 



-3Mr2 + Ma"^ + r^ + ra? 
r A(l + s) - M {a^ -r'^)s 
2 Mrm + As cos 6 

'tp 



ci^2. — ' "-— ' "- , (Al) 

c, ^ -2 '^^^"^^~f ^'^ -' ^' , (A2) 

C3^2a ^"""";;^^^^^\ (A3) 



r -|- ci 
C4 = -2 a m — —: — , (A4) 

(A5) 



the coefficients of M can be written as 



Defining 



db 
nisi = -bci + b —^ + C2 (A6a) 

m32 = 6 C3 - C4 . (A6b) 



A(— ?Ti^ — 2 cos ^sm — cos^ 6' s^ + sin^^s) ,, , 

L^ sm 6 

(r — M) sma , . , 

C6 = -4^ ^^ > (A8) 



(A9) 



the coefficients of A are given by as well as for investigating asymptotic behavior 



Finally, with 



flgi = C5 , (AlOa) 

032 = -C6 , (AlOb) 

ass = ci , (AlOc) 

034 = -C2 . (AlOd) 

(AlOe) 



(All) 
(A12) 



C7 = 


A 

A 

cot e , 


C8 = 



the only non- vanishing coefficient of the operator matrix L reads 



r\2 O 

^31 - C7 ^ + C8 ^ . (A13) 
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TABLES 

TABLE L The numerically determined power-law exponents /x governing the late-time behav- 
ior of |<1>| oc t^^ measured at r* = 20M, = -k 12. The theoretical value for a = is yU = 7. 



ajM 


;U (m = 0) 


//(m = 1) 





7.70 


7.69 


0.25 


7.73 


7.77 


0.50 


7.75 


7.68 


0.75 


7.81 


7.79 


0.90 


7.68 


7.71 



TABLE IL The numerically extracted QNM frequencies that dominate the late-time part of 
the ringing phase for m = 0, ±1. 



ajM 


u)M {m = 0) 


ujM{m= |1|) 





0.373 + i 0.0885 


0.373 + i 0.0885 


0.25 


0.376 + i 0.0885 


0.392 + i 0.0885 


0.50 


0.383 + i 0.0865 


0.420 + i 0.0860 


0.75 


0.398 + i 0.0835 


0.467 + i 0.0795 


0.90 


0.412 + i 0.0780 


0.517 + i 0.0695 
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FIGURES 



m 1 




600 



FIG. 1. Quasinormal ringing for a = 0.99M, M = 0.5, and m = 2. During 120M < t < 160M, 
the field decays as |<5| oc e~^'^'^^^; the subsequent evolution is dictated by a slower decaying mode 
|$| oc 6"°-°'^°^*. 
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FIG. 2. Snapshots of the evolution of ^ji for m = 2, at 6 = tt/2, in the time interval 
< t < 200M. For clarity, ^r/u is displayed with a a scale factor: a = (10, 10^, 10^, 10^) 
for t = (50, 100, 150, 200) M, respectively. The values used for the angular momentum are 
a/M = 0, 0.5, 0.9, 0.99. The initial data consist of a bell-shaped pulse centered at 75M and 
propagating inwards. Differences in the pulses do not become appreciable until the pulse has been 
reflected at the potential barrier outside the black hole (see t = lOOM snapshot). By t = 200M a 
decrease of ~ 25% in the wavelength can be observed. The maximum amplitude of the pulse also 
decreases with a. 
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FIG. 3. Evolution for a/M = 0.99 at t = 50M (a), t = lOOM (b), t = 150M (c), and 
t = 200M (d). The initial data consist of a bell-shaped pulse centered at r* = 75M, 9 = tt/2 and 
propagating inwards. During the evolution there is a clear transition of the angular distribution to 
-2Yi oc cos^ 9/2. 



